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Abstract 

A novel approach is proposed for high-fidelity modeling of progressive damage and failure in composite 
materials that combines the Floating Node Method (FNM) and the Virtual Crack Closure Technique (VCCT) to 
represent multiple interacting failure mechanisms in a mesh-independent fashion. In this study, the approach is 
applied to the modeling of delamination migration in cross-ply tape laminates. Delamination, matrix cracking, 
and migration are all modeled using fracture mechanics based failure and migration criteria. The methodology 
proposed shows very good qualitative and quantitative agreement with experiments. 

1 Introduction 

Damage in composite materials generally occurs as a combination of different and interacting failure 
mechanisms, e.g. delamination and matrix cracking. Capturing these interactions accurately is essential to 
confidently model and predict progressive damage and failure. Many approaches have recently been proposed 
that explicitly model different failure mechanisms and attempt to capture their interaction [1-3]. However, no 
single strategy has yet surfaced as totally satisfactory. The present approach combines the Floating Node 
Method (FNM) [4] with the Virtual Crack Closure Technique (VCCT) [5,6] to explicitly account for different 
failure mechanisms and their interaction. Delamination, matrix cracking, and migration events are all modeled 
with the same FNM element using fracture-mechanics-based failure and migration criteria. The approach is 
validated using recent experimental results in which a setup capable of isolating a single complete migration 
event was developed [7]. Migration is defined as the transition of a delamination from its current ply interface to 
a new interface via transverse ply cracking. The results obtained offer a simple, yet very challenging, 


benchmark for the validation of numerical and analytical approaches aimed at high-fidelity modeling of 
composite materials, such as the one proposed in the present paper. 

2 Floating Node Method (FNM) 

2.1 Introduction 

The FNM is a recently proposed numerical method, capable of representing multiple evolving discontinuities in 
solids [4]. Compared to other methods used for the same purpose, such as extended Finite Element Method 
(XFEM) or the Phantom Node Method (PNM), it provides the following key advantages. First, FNM does not 
introduce an error on the crack geometry when mapping to natural coordinates. Second, it does not require 
numerical integration over only part of a domain. Third, it is ideally suited for the representation and 
incorporation of multiple and complex networks of weak, strong and cohesive, discontinuities. Fourth, FNM 
yields the same solution as a finite element mesh where the discontinuity is represented explicitly, and finally it 
is conceptually simpler than the PNM or XFEM. 

2.2 Element formulation 

The static equilibrium of a body with volume H under body forces with density f (acting on H) and traction t 
acting on the boundary can be expressed in the weak form as: 

J n e T (v)ff(u)dn = J n v T fdfl + J rn v T tdr fi (1) 

where u is the displacement vector; v is the test function; e is the strain tensor related to u through the 
differential operator relative to Cartesian coordinates L x as e = £ x ( u); and a is the stress tensor related to the 
strains through Hooke’s law as a = De, with D being the constitutive tensor. In the Floating Node Method, each 
real node “i” is characterized by its nodal coordinates x t and associated Degrees of Freedom (DoF) q;. In 
addition, an FNM element contains a suitable number of floating DoF without pre-defined associated nodal 
position vectors. These additional floating nodes are used to represent discontinuities. Their number varies with 
the number and type of discontinuities, weak or strong, modeled within each element. A strong discontinuity is 
defined as a jump in a field quantity (e.g. displacements), while a weak discontinuity is defined as a jump in the 
gradient of a field quantity (e.g. strains). Figure 1 shows an example of such an element with four nodes and 
four additional sets of floating DoF required to represent a strong discontinuity. 

2.2.7 Element formulation without weak/strong discontinuity 

If there is no discontinuity to be modeled by the element (e.g. before failure), the formulation is the same as in 
the standard Finite Element Method (FEM). The vector of nodal coordinates is given as x n . In the case of 
Figure 1, x is given by: 
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Xfi = [x 1 ,x 2 ,x 3 ,x 4 ] 

The Jacobian of the transformation from physical coordinates x to natural coordinates ^ is: 


( 2 ) 


dx dN 

, = df = df Xn 


(3) 


where N is a matrix of shape functions. Assuming an isoparametric formulation, the displacement field u is 
related to the real DoF q through: 


u = Nq 


(4) 


In the case of Figure 1, q T = [q 1; q 2 , ({ 3 , ({ 4 ]. The strains can also be expressed in terms of the real DoF q as: 


e = £ x (u) = ^(N)J *q = Bq 
where B = £^(N)J -1 , leading to the stiffness matrix: 
K = J B t D Bdet(J)dS 


(5) 


( 6 ) 


where E is the integration domain (in natural coordinates) corresponding to kl (in physical coordinates), and to 
the vector of nodal forces: 


Q = J N T fdet(J)d E + J N T fdet(J)dS 


(7) 


where is similarly the boundary corresponding to r n . The weak form of equilibrium, Equation 1, becomes: 

Kq = Q (8) 

2.2.2 Element formulation with weak/strong discontinuity 

Once a discontinuity in the element is defined, the element is split in two or more partitions (depending on the 
discontinuity). Without loss of generality, a case in which the element is split in two partitions, kl A and H B , is 
illustrated in Figure 1. For each partition, kl A and H B , a vector of nodal coordinates, x^ andx^ B , is defined. 
For the case in Figure 1 : 


X L = [x 5 ,X 6 ,X 3 ,X 4 ] 


(9) 


xL = [x 1 ,x 2 ,x 6 ,x 5 ] 


( 10 ) 
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Each partition has a separate Jacobian: 


dx dN 


( 11 ) 


,A d? d? XflA 


dx dN 


( 12 ) 


,B d? d? XnB 


The displacements u A and u B , in partitions H A and H B , respectively, are interpolated separately from the DoFs 
q A and q B associated with each partition: 


These DoFs, q A and q B , did not have an associated position beforehand, and therefore were considered to be 
“floating”. As the discontinuity is defined, they are linked to the respective position vector. In the case of Figure 
1, which represents a strong discontinuity, qj = [q 5 , q 6 , q 3 , q 4 ] and qj = [q 1; q 2 , qs,] . Note that there are 
four different sets of floating DoF: q 5 is different from q 5 , and q 6 is different from q 6 ,. If a weak discontinuity 
were to be modeled, only two sets of floating DoF would be included in the element, and the DoF with a prime 
would coincide with those without a prime. 

The strains then become: 


u A = Nq A and u A = Nq A 


( 13 ) 


e A = £ x (u a ) = ^(N)J a *q A = B A q A 


( 14 ) 


e B — ^x( u b) — £^(N)Jb — ^bQb 


( 15 ) 


The stiffness matrices for partitions fl A and fl B are: 



( 16 ) 



( 17 ) 


and the force vectors: 


Q A = 4 N T f det(J A )d2 + ; rs N T f det(I A )dS 
Q b = 4 N T f det(J B )dS + 4 b N T f detO B )dS 


( 18 ) 


( 19 ) 
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Finally, the equations of equilibrium can be written as: 

KaQa = Qa and K B q B = Q B 


2.2.3 Element topology and assembly 

To illustrate this approach, a floating node element capable of modeling two weak interfaces with an arbitrary 
crack between them was implemented, as shown in Figure 2(a). The nodal position of the initial integration 
domain, and the real and floating DoF are indicated in the figure. Each floating DoF is associated with either an 
edge or the inner domain of the element. Once weak/strong discontinuities are detected, the floating DoF at the 
edges are used to determine the solution for the nodal positions given by the intersection of these discontinuities 
with the element edges (Figs. 2(b) to 2(e)). The additional inner floating DoF are used to determine the solution 
for the nodal positions originated by the intersection of multiple cracks within the element (Fig. 2(f)). In Figures 
2(b) to (f) only the floating DoF used are represented. In each case, the floating DoF that are not used have no 
assigned position, only a topological relation to an edge or to the inner domain. During assembly, inner floating 
DoF can be removed from the analysis through static condensation. The floating DoF associated with the edges 
are assembled with the corresponding DoF of neighboring elements, and will therefore have a unique position in 
the global DoF vector. All floating DoF are then used as required to model weak/strong discontinuities as they 
evolve throughout the simulation. 

The approach can be extended to model an arbitrary number of weak and strong discontinuities within an 
element, and therefore an arbitrary number of plies and interacting cracks, provided the topology is adequately 
defined, and a sufficient number of floating DoF are used. 

3 Virtual Crack Closure Technique (VCCT) 

In the present work, the FNM method has been coupled to the Virtual Crack Closure Technique (VCCT) [5]. An 
extensive review of VCCT is provided in [6]. Traditionally, VCCT is used to obtain energy release rates in cases 
where the crack path is known beforehand and can be aligned with the elements’ edges, such as the case for 
delamination. 

In VCCT, the Mode I and II energy release rates are obtained from [6]: 
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where F n and F t are the normal and tangential components of internal force vector F at the crack tip; lq n 1 and 
lq t J are the normal and tangential components of the displacement jump [q] between the nodes immediately 
behind the crack tip; and a x and a 2 are the lengths of the crack in the elements behind and ahead of the crack 
tip, respectively (Fig. 3). When combining VCCT with FNM, the forces and displacements needed to obtain the 
energy release rates, Eqs (21) and (22), are computed at the floating DoF as the crack develops (Fig. 3). 
Additionally, floating DoF are also added along a virtual crack plane up to a distance r, named the ‘enrichment 
radius’, r enr . The associated floating nodes provide additional degrees of freedom along this path, which enable 
the accurate modeling of the deformations near the crack tip [4]. At the end of the enrichment radius, the 
displacements of the floating DoF are interpolated from the displacements of the real DoF. In the present work, 
the enrichment radius was chosen to be equal to the largest in-plane dimensions of the numerical model. 

4. Delammation propagation 

In this paper, delamination is modeled using VCCT to compute the energy release rates Gj and G n at each 
delamination front position. These are then used in a failure criterion: 

fiGM =^F-1 = 0 (23) 

G C 

where G T = Gj + G Ih and G c is the critical energy release rate given by [8]: 

(GuX 71 

G c ~ G IC + (fine ~ G/c) \^~) (24) 

For delamination growth according to this criterion, G c is assumed to equal the critical energy release rate of the 

interface, Gc nt . The delamination front is advanced by one element along the projected crack path when 

r — r'/nt 
Ur — Li c . 

5 Matrix crack propagation 

In the present study, matrix crack propagation is also modeled with VCCT. However, unlike delamination, the 
matrix crack is assumed to propagate following a Mode I path in the through- thickness direction, as supported 
by experimental evidence [7]. Therefore, the failure criterion is assumed to be given by: 

/(G/) = 1 = 0 (25) 

Gc 

where the critical energy release rate is given by G c = G lc . 

At each crack growth increment, the Mode I crack path orientation is approximately determined using a 
maximum tangential stress criterion. This criterion, typically written in terms of the Mode I and II stress 
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intensity factors [9], can also be written in terms of energy release rates Gj and G Ih where the angle 6 that 
maximizes the tangential stress Gqq is obtained by evaluating: 



and for the two solutions obtained choosing the one that maximizes Oqq given by: 


oee = Eh (3 cos (0 

+ cos (t)) + 

./G^Sgn(F t ) [-3 sin 

(9 -3 *»(¥)} 


(26) 


(27) 


where F t is the tangential component of the internal force vector F at the crack tip. The tangential stress Oqq 
relates to oqq by: 




where, 



plane stress 
plane strain 


(29) 


(30) 


E is the Young’s modulus; and v the Poisson’s ratio. Once the angle is determined, the criterion given by Eq. 
(25) is assessed. If the criterion is met, the crack advances by one element along the projected crack path. 

6 Migration 


6.1 Delamination to matrix crack - migration criterion 


Extensive research has been done on migration of cracks from a bi-material interface. He and Hutchinson [10] 
have proposed a necessary condition for migration to occur: 


r , ts r 

u max ^ U T 

„ D 


Int 


Int 


(31) 


This condition states that for a crack to migrate from an interface into a given Material B, as shown in Figure 4, 
the ratio between the energy release rate associated with its growth into that material G^ ax , (at an angle co max 
that maximizes its value) and the critical energy release rate of that material , has to be higher than the ratio 
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between the energy release rate associated with the growth of the crack along the interface Gf nt , and the critical 
energy release rate G c nt of that interface. 

With the approach proposed, evaluating Eq. (31) would require more than one FE analysis at each interface 
crack position. In each of these analyses, a branched crack would need to be initiated at different tentative 
angles co in order to determine the angle, 0 ) max , that would maximize G B . 

In the present work, an alternative criterion is proposed. With this criterion, all variables needed to determine 
whether migration occurs, can be directly obtained from the FNM results at the each interface crack position. 
The criterion states that: (i) the sign of the tangential component of the internal force vector F t , associated with 
Mode II shear displacement, dictates the direction and therefore the material into which the interface crack 
migrates and (ii) migration can only be completed if it is energetically favorable. The criterion is given by 
1, 0,-1} as: 


Mr = 


c (G' T nt G!r nt \ , 

Sgn h?--c?f) +1 


Sgll(F t ), Sgll(F t ) < 0 
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, (C,f‘ c,';“ i , ,\ 

Sgn (-er-c=J +1 l 


Sgn(F t ), Sgn(F t ) > 0 


/ 


(32) 


If Mb = 0 the interface crack does not migrate, and if = 1 or = — 1 the crack migrates into material A 
or B, respectively. In the present form, Eq. (32) assumes that internal force vector is computed at the lower 
surface (Material B). If the internal force vector is computed at the upper surface the signs of the inequalities in 
Eq. (32) need to be reversed. Depending on the sign of F t , the ratio between the energy release rate determined 

Qlnt Q^ nt 

at the interface, and the critical energy release rate of Material A or B,-^-or-^-, is compared to the ratio 
between the energy release rate determined at the interface and the critical energy release rate of the interface, 

Qlnt G^ nt G^ nt G^ nt 

If the ratio or is greater than the necessary condition for migration is met, and = 1 or 
g c g c G c G c 

= —1. With this criterion, all variables needed to determine whether migration occurs, F t and Gf nt , can be 
directly obtained from the FNM results. Note that although this criterion is necessary, it is not sufficient for 
migration to occur unless the ratio being assessed is greater than one. 

In the present work, the interfacial crack, or delamination, propagates between a 0° and a 90° ply. Assuming 
Material A is the 0° ply, and Material B is the 90° ply, the critical energy release rate to propagate a crack into a 
0° ply is much greater than the critical energy release rate of the interface, since fiber fracture is required for that 
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o G^ nt G^ nt 

process to occur, i.e. G<?° — » &c nt [H]- Thus, is always smaller than Therefore, even if the sign of 

U Q G Q 

F t favors migration into Material A, migration will not occur, since it is not energetically favorable. Typically 
for composite laminates, the critical energy release rate of an interface crack (delamination) and of a matrix 
crack are of the same order of magnitude, i.e. G<?°° ~ 6 c nt [12]. Hence, Gc nt is assumed to be given by Eq. (24) 
and Gf °° to be equal toG /c . Therefore, assuming Material B is the 90° ply, provided the sign of F t favors 
migration into Material B and it is energetically favorable, migration can occur. 

6.2 Matrix crack to delamination - migration criterion 

When a matrix crack reaches a weak interface, it is assumed to trigger delamination. This delamination is at first 
contained within the FNM element, as illustrated in Figure 2(f). In the following steps, the delamination will 
propagate, or migrate, as detailed in Sections 4 and 6.1, respectively. 

7 Simulation of delamination migration specimens 

7.1 Delamination migration test configuration 

In reference [7] a test aimed at investigating delamination migration was proposed. Its configuration is 
illustrated in Figure 5. The specimen possesses three main features that permit the controlled observation of 
delamination growth followed by migration to another ply interface. First, the specimen geometry is in the form 
of a beam with the intent of promoting uniform delamination growth and migration across the specimen width. 
Second, the specimen contains a PTFE insert (acting as an artificial delamination) at an interface between a 0° 
ply (specimen span direction) and a stack of four 90-degree plies (specimen width direction) (Fig. 5). This 
provides an opportunity for the delamination to migrate to another ply interface by kinking through the 90° 
degree ply stack. Third, the specimen can be loaded in a manner to cause delamination growth from the PTFE 
insert that eventually migrates to another ply interface. This sequence of fracture events is made possible by the 
way in which specimen loading affects shear stresses acting across the delamination front. 

7.2 Finite Element Model 

The experiments presented in [7] were simulated using the finite element solver ABAQUS/Standard® (Implicit). 
Plane-stress and a small-displacement formulation were used. The floating node element (Fig. 2) was developed 
and implemented as a user defined element (UEF) in ABAQUS, and applied in the center region of the model, 
as shown in Figure 6. The specimen layup used in the experiments and in the model, is also given in Figure 6. In 
the layup description T represents the PTFE insert. Each block of plies of the same orientation was modeled 
with a separate element CPS4 (through-thickness), except for the highlighted plies at the center of the model 
(Fig. 6). These were modeled within a single FNM element, as illustrated in Figure 2. The material properties 
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used are given in Tables 1 and 2. A prescribed displacement was applied at different positions, L, at a single 
node, to simulate the displacement-controlled test. Four load positions were simulated corresponding to 
L = { a 0 , l.la 0 , 1.2a 0 , 1.3a 0 } (Fig. 6). In the experiments, the specimens were subject to moisture intake from 
the ambient air leading to relaxation of residual thermal stress. Therefore, residual thermal stresses were not 
considered. 

7.3 Load-displacement response and sequence of events 

The computed load-displacement curves obtained for different values of L are given in Figure 7. Even 
accounting for the compliance of the test fixture, measured in [7], all load-displacement curves obtained from 
the simulations show a slightly stiffer response than what was observed experimentally. The maximum load is 
in general very well predicted for the cases L = a 0 to L = 1.2a 0 , and slightly over-predicted for the case 
L = 1.3 a 0 (Figs. 7(a) to 7(d)). For L = a 0 , the simulations predict unstable delamination after the peak load 
(Fig. 7(a)). This unstable delamination stops just before migration. After the migration event, delamination 
continues stably along the next 0°/90° interface. A similar sequence of events was observed experimentally [7]. 
However, in the experiments, after the first unstable event, further loading was necessary before migration 
occurred. For the cases L = 1.1 a 0 to L = 1.3 a 0 , the simulations predict that (i) the first peak in load is followed 
by a region of stable delamination growth (Figs 7(b) to 7(d)) that increases with the load offset, and that (ii) 
stable delamination is followed by a region of unstable delamination, which corresponds to the load-drop in 
Figures 7(b) to 7(d). The migration event occurs within this load-drop. After the migration, unstable 
delamination continues propagating along the next 0°/90° interface. Loading the specimen further eventually 
leads to a transition from unstable to stable delamination propagation. This sequence of events is very similar to 
what was reported experimentally [7]. The main difference is that in the experiments, for the cases L = 1.2 a 0 
andL = 1.3a 0 , after the peak load, delamination propagates through a series of unstable events, rather than 
stably, as predicted by the simulations. 

7.4 Migration location 

Figure 8(a) provides an illustration of the typical crack paths predicted by the simulations and obtained 
experimentally. Figure 8(b) compares numerical and experimental results for the distance from the load- 
application point at which migration onset occurs, A K , as a function of the normalized load offset, L/a 0 . 
Overall, except for L = a 0 , migration is predicted to occur closer to the load application point than what was 
observed experimentally. The simulations also predict a decreasing trend in A K with L/a 0 , which is similar to 
what was observed experimentally for the cases L = 1.1 a 0 to L = 1.3a 0 . 
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7.5 Matrix crack path 

In Figure 9, the simulated matrix crack paths obtained for the different load positions was superimposed on an 
image of a typical crack path observed experimentally, where L = a 0 . The predicted matrix crack path agrees 
well the crack path observed experimentally [7]. The simulations do not predict variation in the matrix crack 
path with a change in load offset. 

7.6. Discussion 

All load-displacement curves obtained numerically show a slightly stiffer response than the experimental 
results. In the simulations, the engineering properties used were obtained from uniaxial test data. These are 
known to be higher than their flexural counterparts [13]. Since the specimen is predominantly loaded in flexure, 
the difference between uniaxial and flexural properties is therefore likely to be the cause for the overall stiffer 
response obtained in the simulations. Overall, the simulated load-displacement curves, sequence of events, and 
failure morphology are in very good agreement with the experimental results. For the case L = a 0 , 
experimental results show an increase in load at the end of the first unstable event, before migration. This 
increase in load is likely to be caused by the development of fiber bridging during the unstable part of the 
delamination, which was not accounted for in the model. For the load cases L = 1.2 a 0 to L = 1.3 a 0 , the only 
significant difference between simulations and experiments is that delamination growth onset and propagation, 
prior to the final load-drop see (Figs. 7(c) and 7(d)), does not occur stably as predicted in the simulations, but 
rather through a series of unstable events. In these experiments the unstable delamination growth onset was 
likely due to the presence of resin rich pockets at the tip of the PTFE insert, while the subsequent unstable 
events, prior to the final load-drop, may be due to features accompanying delamination growth, such as fiber 
bridging. Both of these were not accounted for in the simulations. 

The migration location was slightly under-predicted for all cases except L = a 0 . Nevertheless, the migration 
location trend observed experimentally, for the cases L = 1.1 a 0 to L = 1.3 a 0 , was well captured. This further 
supports the validity of the approach, and in particular, of the migration criterion proposed. The latter is 
additionally confirmed by the good agreement between the predicted matrix crack path and that observed 
experimentally. This also suggests that the criterion used to determine the matrix crack orientation (as it 
propagates in the through-thickness direction) is adequate. 
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8. Conclusions 


A novel approach that combines the FNM (Floating Node Method) with VCCT (Virtual Crack Closure 
Technique) has been demonstrated to adequately capture delamination migration in cross-ply laminates. 
Delamination, matrix cracking, and migration, are all modeled with the same FNM element using fracture 
mechanics based failure and migration criteria. Results show very good agreement with the maximum load, and 
migration location obtained experimentally. In addition, simulations capture the overall failure morphology and 
typical crack path, observed in the experiments. 
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Table 1. Elastic properties, IM7-8552 [14]. 


Fn 

^22 — £*33 

V 12 = v 13 

v 23 

6l2 - 6l3 

^23 

a ± 

a 2 = a 3 

161.0 (GPa) 

11.38 (GPa) 

0.32 

0.436 

5.17(GPa) 

3.98(GPa) 

lxlO -5 K -1 

31xl0 -6 K -1 


Table 2. Critical energy release rates for delamination of IM7-8552 [15]. 


G/c 

Giic V 

0.21 (kj/m 2 ) 

0.77(kJ/m 2 ) 2.1 


Without discontinuity With weak/strong discontinuity 





Key: 
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Figure 1. Overview of the Floating Node Method (FNM). 
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(a) topology of the floating node element 



Key: 

O Nodal coordinates of 
initial integration 
domain 

□ Nodal coordinates of 
crack position 

a Degree of freedom 

A Floating degree of 
freedom 






crack 


(f) matrix crack with 
migration to delamination 


Figure 2. Floating node element used showing different weak/strong discontinuities scenarios 




Figure 3. Virtual Crack Closure Technique and 
Floating Node Method for arbitrary crack 
propagation. 


Figure 4. Bi-material interface. 






